//==============================================================================
//
//           DComplex - the complex number class in the ofc-library
//
//               Copyright (C) 2003  Dick van Oudheusden
//  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public
// License along with this library; if not, write to the Free
// Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//
//==============================================================================
// 
//  $Date: 2005-03-29 18:22:58 $ $Revision: 1.11 $
//
//==============================================================================

#include "ofc/DComplex.h"

#include <math.h>

#if _INTERFACE_

#include "ofc/config.h"

#include "ofc/DTextable.h"
#include "ofc/DSimpleArithmetic.h"

#include "ofc/DText.h"

// The complex class implements a number of methods for working with complex numbers.
//
// Heavily 'inspired' by the excellent ccmath library by Daniel A. Atkinson (c) Copyright 2000.

@interface DComplex : Object <DTextable,DSimpleArithmetic>
{
@private
  double       _re;
  double       _im;
}

#endif



@implementation DComplex


//// Constructors

//
// Initialise to complex number zero
//
// @return the object
//

- (DComplex *) init
{
  [super init];

  _re = 0.0;
  _im = 0.0;
  
  return self;
}

//
// Initialise to a complex number
//
// @param re       the real part
// @param im       the imaginary part
//
// @return the object
//

- (DComplex *) init :(double) re :(double) im
{
  [self init];

  _re = re;
  _im = im;

  return self;
}



//// Member methods

//
// Get the real part
//
// @return the real part
//

- (double) re
{
  return _re;
}

//
// Set the real part
//
// @param re       the real part of a complex number
//
// @return the object
//

- (DComplex *) re :(double) re
{
  _re = re;

  return self;
}

//
// Get the imaginary part of the number
//
// @return the imaginary part
//

- (double) im
{
  return _im;
}

//
// Set the imaginary part of the number
//
// @param im       the imaginary part
//
// @return the object
//

- (DComplex *) im :(double) im
{
  _im = im;

  return self;
}

//
// Set the real and imaginary part of the number
//
// @param re       the real part
// @param im       the imaginary part
//
// @return the object
//

- (DComplex *) set :(double) re :(double) im
{
  _re = re;
  _im = im;

  return self;
}

//
//
// Move an other object in this object
//
// @param other    the other object
//
// @return this object
//

- (DComplex *) move :(DComplex *) other
{
  _re = [other re];
  _im = [other im];

  return self;
}



//// Methods on complex numbers

//
// Add with a complex number 
//
// @param other    the complex number
//
// @return the object
//

- (DComplex *) add :(DComplex *) other
{
  _im += [other im];
  _re += [other re];
  
  return self;
}

//
// Add two complex numbers and store the result in the object
//
// @param src1     the first complex number
// @param src2     the second complex number
//
// @return the object
//
- (DComplex *) add :(DComplex *) src1 :(DComplex *) src2
{
  _im = [src1 im] + [src2 im];
  _re = [src1 re] + [src2 re];
  
  return self;
}

//
// Subract with a complex number
//
// @param other     the complex number
//
// @return the object
//

- (DComplex *) sub :(DComplex *) other
{
  _im -= [other im];
  _re -= [other re];
  
  return self;
}

//
// Subract two complex numbers and store the result in the object
//
// @param src1      the first complex number
// @param src2      the second complex number
//
// @return the object
//
- (DComplex *) sub :(DComplex *) src1 :(DComplex *) src2
{
  _im = [src1 im] - [src2 im];
  _re = [src1 re] - [src2 re];
  
  return self;
}

//
// Multiply with an other complex number
//
// @param other    the complex number
//
// @return the object
//

- (DComplex *) mul :(DComplex *) other
{ 
  double re = [other re];
  double im = [other im];
  double tm;

   tm = (_re * re - _im * im);
  _im = (_re * im + re * _im);
  _re = tm;

  return self;
}

//
// Multiply two complex numbers and store the result in the object
//
// @param src1     the first complex number
// @param src2     the second complex number
//
// @return the object
//
- (DComplex *) mul :(DComplex *) src1 :(DComplex *) src2
{ 
  double re1 = [src1 re];
  double im1 = [src1 im];
  double re2 = [src2 re];
  double im2 = [src2 im];

  _re = (re1 * re2 - im1 * im2);
  _im = (re1 * im2 + re2 * im1);

  return self;
}

//
// Multiply with a real number
//
// @param re       the real number
//
// @return the object
//

- (DComplex *) rmul :(double) re
{
  _re *= re;
  _im *= re;

  return self;
}

//
// Multiply with a imaginary number
//
// @param im       the imaginary number
//
// @return the object
//

- (DComplex *) imul :(double) im
{
  double tm;

   tm =  _re * im;
  _re = -_im * im;
  _im =   tm;

  return self;
}

//
// Divide with a complex number
//
// @param other    the divider
//
// @return the object
//
- (DComplex *) div :(DComplex *) other
{
  double re = [other re];
  double im = [other im];
  double tm;
  double dv;

   dv =   re * re +  im * im;
   tm = (_re * re + _im * im) / dv;
  _im = (_im * re - _re * im) / dv; 
  _re = tm;

  return self; 
}

//
// Divide two complex numbers and store the result in the object
//
// @param src1     the first complex number
// @param src2     the divider
//
// @return the object
//
- (DComplex *) div :(DComplex *) src1 :(DComplex *) src2
{
  double re1 = [src1 re];
  double im1 = [src1 im];
  double re2 = [src2 re];
  double im2 = [src2 im];
  double dv;

   dv =  re2 * re2 + im2 * im2;
  _re = (re1 * re2 + im1 * im2) / dv;
  _im = (im1 * re2 - re1 * im2) / dv; 

  return self; 
}

//
// Conjugate the complex number
//
// @return the object
//

- (DComplex *) cng 
{
  _im = -_im;

  return self;
}

//
// Calculate the modulus of the complex number
//
// @return the modulus
//

- (double) abs 
{
  return sqrt(_re * _re + _im * _im);
}

//
// Calculate the square of the modulus of the complex number
//
// @return the square of the modulus
//

- (double) nrm
{
  return (_re * _re + _im * _im);
}

//
// Calculate the square root for the complex number
//
// @return the object
//

- (DComplex *) sqrt
{
  double r = [self abs];

  r = sqrt(ldexp(r + fabs(_re), -1));
  if (r == 0.0) 
  {
    _re = 0.0;
    _im = 0.0;
  }
  else
  {
    if (_re >= 0.0)
    { 
      _re = r; 
      _im = ldexp(_im / r, -1);
    }
    else
    { 
      _re = ldexp(fabs(_im) / r, -1);
      if (_im >= 0.0) 
      {
        _im = r; 
      }
      else
      {
        _im= -r;
      }
    }
  }
  return self;
}

//
// Calculate the exponent function for the complex number
//
// @return the object
//

- (DComplex *) exp
{
  double s = sin(_im);
  double c = cos(_im);
  double u = exp(_re);

  _re = u * c; 
  _im = u * s;

  return self;
}

//
// Calculate the natural logarithm for the complex number
//
// @return the object
//

- (DComplex *) log
{
  double r = [self nrm];

  _im = atan2(_im, _re);
  _re = ldexp(log(r), -1);

  return self;
}

//
// Calculate the trigonometric functions sine for the complex number
//
// @return the object
//

- (DComplex *) sin
{
  double s = sin(_re);
  double c = cos(_re);
  double u = exp(_im); 
  double v = 1.0 / u;

  u = ldexp(u + v, -1); 
  v = u - v;

  _re = u * s;
  _im = c * v;

  return self;
}

//
// Calculate the trigonometric functions cosine for the complex number
//
// @return the object
//

- (DComplex *) cos
{
  double s = sin(_re);
  double c = cos(_re);
  double u = exp(_im); 
  double v = 1.0/u;

  u = ldexp(u + v, -1);
  v = u - v;

  _re =  c * u;
  _im = -s * v;

  return self;
}

//
// Calculate the trigonometric functions trangent for the complex number
//
// @return the object
//

- (DComplex *) tan
{
  double s = sin(_re);
  double c = cos(_re);
  double u = exp(_im); 
  double v = 1.0 / u;
  double d;

  u = ldexp(u + v, -1);
  v = u - v;
  d = c * c + v * v; 

  _re = s * c / d; 
  _im = u * v / d;

  return self;
}

//
// Calculate the inverse trigonometric function sine for the complex number
//
// @return the object
//

- (DComplex *) asin
{
  DComplex *u = [DComplex alloc]; 
  DComplex *w = [self copy];
  
  double    t;

  [u init :1.0 :0.0];
  [w move :self];
  
  [w mul :w];
  [u sub :w];
  [u sqrt  ];

  t   = _re;  
  _re = [u re] - _im; 
  _im = [u im] +  t;
  
  [self log];

  t   = _re;
  _re = _im; 
  _im = -t;

  [u free];
  [w free];
  
  return self;
}

//
// Calculate the inverse trigonometric function cosine for the complex number
//
// @return the object
//

- (DComplex *) acos
{
  DComplex *u = [DComplex alloc];
  DComplex *w = [self copy];
  double    t;

  [u init :1.0 :0.0];
  
  [w mul :w];
  [u sub :w];
  [u sqrt  ];
  
  _re -= [u im]; 
  _im += [u re]; 

  [self log];

   t  = _re;
  _re = _im; 
  _im =  t;

  return self;
}

//
// Calculate the inverse trigonometric function tangent for the complex number
//
// @return the object
//

- (DComplex *) atan
{
  DComplex *u = [DComplex alloc];
  DComplex *w = nil;

  [u init :-_im :_re];
  
  _re = 1.0;
  _im = 0.0;

  w = [self copy];

  [self add :u];
  [w    sub :u];
  [self div :w];

  [self log];
  [self imul :-0.5];

  return self;
}

//
// Calculate the hyperbolic function sine for the complex number
//
// @return the object
//

- (DComplex *) sinh
{
  double s = sin(_im);
  double c = cos(_im);
  double u = exp(_re);
  double v = 1.0 / u;

  u = ldexp(u + v, -1);
  v = u - v;

  _re = v * c;
  _im = u * s;

  return self;
}

//
// Calculate the hyperbolic function cosine for the complex number
//
// @return the object
//

- (DComplex *) cosh
{
  double s = sin(_im);
  double c = cos(_im);
  double u = exp(_re);
  double v = 1.0 / u;

  u = ldexp(u + v, -1);
  v = u - v;

  _re = c * u;
  _im = v * s;

  return self;
}

//
// Calculate the hyperbolic function tangent for the complex number
//
// @return the object
//

- (DComplex *) tanh
{
  double s = sin(_im);
  double c = cos(_im);
  double u = exp(_re); 
  double v = 1.0 / u;
  double d;

  u = ldexp(u + v, -1); 
  v = u - v;
  d = c * c + v * v; 

  _re = u * v / d;
  _im = s * c / d;

  return self;
}

//
// Calculate the inverse hyperbolic function sine for the complex number
//
// @return the object
//

- (DComplex *) asinh
{ 
  DComplex *u = [DComplex alloc];
  DComplex *w = [self copy];

  [u init :1.0 :0.0];
  
  [w mul :w];
  [u add :w];
  [u sqrt  ];

  [self add :u];
  [self log   ];

  [u free];
  [w free];

  return self;
}

//
// Calculate the inverse hyperbolic function cosine for the complex number
//
// @return the object
//

- (DComplex *) acosh
{
  DComplex *u = [DComplex alloc];
  DComplex *w = [self copy];
  int      kf = ((_im == 0) && (_re < -1)) ? 1 : 0;

  [u init :1.0 :0.0];
  
  [w mul  :w];
  [w sub  :u];
  [w sqrt   ];

  [u move :w];

  [self add :u];
  [self log   ]; 

  if (_re < 0.0)
  {
    _re= -_re;
    _im= -_im;
  }
  if (kf) 
    _im= -_im;

  [u free];
  [w free];

  return self;
}

//
// Calculate the inverse hyperbolic function tangent for the complex number
//
// @return the object
//

- (DComplex *) atanh
{
  DComplex *u = [DComplex alloc]; 
  DComplex *w = [DComplex alloc]; 

  [[w init :1.0 :0.0] sub :self];
  [[u init :1.0 :0.0] add :self];

  [u div :w   ];
  [u log      ];
  [u rmul :0.5];

  [self move :u];

  [u free];
  [w free];

  return self;
}


//// Textable protocol methods

//
// Convert to a text string 
//
// @return a (new) text string with the complex number
//

- (DText *) toText
{
  DText *str = [[DText alloc] init];

  if (_re == 0.0)
  {
    [str format :"%gj",_im];
  }
  else if (_im == 0.0)
  {
    [str format :"%g", _re];
  }
  else if (_im < 0.0)
  {
    [str format :"%g%gj",_re,_im];
  }
  else
  {
    [str format :"%g+%gj",_re,_im];
  }
  
  return str;
}

@end

/*===========================================================================*/

